# Cell 1: Imports and configuration
import os
import math
import cv2
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.cluster import DBSCAN
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_error

# User settings
TOTAL_INPUT_PATH = r"C:\Users\adith\Downloads\2mM TBPB_.jpg"  # <- change if needed
INNER_SCALE = 0.72
EXPECTED_COLS = 12
CONCENTRATIONS = [0.5,1,2,3,4,5,6,7,8,9,10]
SATURATION_THRESHOLD = 35
MIN_BLOB_AREA = 60
PLOT_DPI = 100

# Create output folder for debugging images if you want to save (optional)
OUT_DEBUG = None  # set to a path string to save images, e.g. r"D:\FYPROJJ\debug"
if OUT_DEBUG:
    os.makedirs(OUT_DEBUG, exist_ok=True)

# Display defaults
plt.rcParams["figure.dpi"] = PLOT_DPI

# Cell 2: Utility functions

def imshow_rgb(img_bgr, title=None, figsize=(8,6)):
    """Show BGR image in notebook as RGB with matplotlib"""
    plt.figure(figsize=figsize)
    plt.imshow(cv2.cvtColor(img_bgr, cv2.COLOR_BGR2RGB))
    if title:
        plt.title(title)
    plt.axis('off')
    plt.show()

def imshow_gray(img_gray, title=None, cmap='gray', figsize=(8,6)):
    plt.figure(figsize=figsize)
    plt.imshow(img_gray, cmap=cmap)
    if title:
        plt.title(title)
    plt.axis('off')
    plt.show()

def to_hsv(img_bgr):
    return cv2.cvtColor(img_bgr, cv2.COLOR_BGR2HSV)

def hsv_to_display_rgb(hsv):
    # convert HSV back to RGB for pleasant display
    return cv2.cvtColor(hsv, cv2.COLOR_HSV2RGB)

def mask_from_hsv(hsv, s_thresh=30, v_thresh=40):
    s = hsv[:,:,1]
    v = hsv[:,:,2]
    mask = ((s > s_thresh) & (v > v_thresh)).astype(np.uint8) * 255
    return mask

def morphological_clean(mask, open_k=3, close_k=5):
    ko = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (open_k,open_k))
    kc = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (close_k,close_k))
    m_open = cv2.morphologyEx(mask, cv2.MORPH_OPEN, ko, iterations=1)
    m_close = cv2.morphologyEx(m_open, cv2.MORPH_CLOSE, kc, iterations=1)
    return m_open, m_close

def contours_from_mask(mask):
    cnts, _ = cv2.findContours(mask.copy(), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    return cnts

def contours_overlay(img_bgr, contours, color=(0,255,0), thickness=2, show=True, title=None):
    out = img_bgr.copy()
    cv2.drawContours(out, contours, -1, color, thickness)
    if show:
        imshow_rgb(out, title=title)
    return out

def contours_to_circles(contours, area_min=MIN_BLOB_AREA, circ_min=0.3):
    """Return list of (x,y,r,area,circularity) for contours passing thresholds."""
    blobs = []
    for c in contours:
        area = cv2.contourArea(c)
        if area < 5:
            continue
        peri = cv2.arcLength(c, True)
        if peri == 0:
            continue
        circ = 4*math.pi*area/(peri*peri)
        (x,y), r = cv2.minEnclosingCircle(c)
        if area >= area_min and circ >= circ_min:
            blobs.append((int(x), int(y), int(round(r)), float(area), float(circ)))
    # left-to-right
    blobs = sorted(blobs, key=lambda b: b[0])
    return blobs

def overlay_circles(img_bgr, blobs, inner_scale=INNER_SCALE, annotate=True, title=None):
    out = img_bgr.copy()
    for i,(x,y,r,area,circ) in enumerate(blobs):
        cv2.circle(out, (x,y), int(r), (0,255,0), 2)              # outer
        cv2.circle(out, (x,y), max(1,int(r*inner_scale)), (255,0,0), 2)  # inner ROI
        if annotate:
            cv2.putText(out, str(i+1), (x-10,y+int(r)+12), cv2.FONT_HERSHEY_SIMPLEX, 0.5, (0,0,0), 2)
    imshow_rgb(out, title=title)
    return out

def hough_circles_fallback(img_bgr, dp=1.2, minDist=24, param1=80, param2=28, minR=12, maxR=60):
    gray = cv2.cvtColor(img_bgr, cv2.COLOR_BGR2GRAY)
    blur = cv2.GaussianBlur(gray, (5,5), 1.2)
    circles = cv2.HoughCircles(blur, cv2.HOUGH_GRADIENT, dp=dp, minDist=minDist,
                               param1=param1, param2=param2, minRadius=minR, maxRadius=maxR)
    if circles is None:
        return []
    circles = np.round(circles[0]).astype(int)
    blobs = [(int(x),int(y),int(r),0.0,0.0) for (x,y,r) in circles]
    blobs = sorted(blobs, key=lambda b: b[0])
    return blobs

def extract_inner_rgb_and_saturation(img_bgr, img_hsv, x, y, r, inner_scale=INNER_SCALE):
    inner_r = max(1, int(r * inner_scale))
    h, w = img_bgr.shape[:2]
    mask = np.zeros((h,w), dtype=np.uint8)
    cv2.circle(mask, (x,y), inner_r, 255, -1)
    px = img_bgr[mask==255]
    hsv_px = img_hsv[mask==255]
    if px.size == 0 or hsv_px.size == 0:
        return (np.nan, np.nan, np.nan, np.nan, np.nan)
    B = float(np.mean(px[:,0])); G = float(np.mean(px[:,1])); R = float(np.mean(px[:,2]))
    S_mean = float(np.mean(hsv_px[:,1])); V_mean = float(np.mean(hsv_px[:,2]))
    RGB_mean = (R+G+B)/3.0
    return R, G, B, RGB_mean, S_mean

def cluster_rows_by_y_from_blobs(blobs, eps_mul=1.8, min_samples=2):
    """Cluster by Y coordinate and return ordered list of lists of blob tuples."""
    if len(blobs) == 0:
        return []
    ys = np.array([b[1] for b in blobs]).reshape(-1,1).astype(float)
    radii = np.array([b[2] for b in blobs])
    med_r = max(6.0, np.median(radii))
    eps = med_r * eps_mul
    db = DBSCAN(eps=eps, min_samples=min_samples).fit(ys)
    labels = db.labels_
    label_map = {}
    for lbl,b in zip(labels, blobs):
        if lbl == -1:
            continue
        label_map.setdefault(int(lbl), []).append(b)
    rows = []
    for lbl,blist in label_map.items():
        avg_y = np.mean([b[1] for b in blist])
        rows.append((lbl, avg_y, blist))
    rows_sorted = sorted(rows, key=lambda t: t[1])
    ordered = [r[2] for r in rows_sorted]
    return ordered

# Cell 3: Full-image stepwise pipeline outputs

img = cv2.imread(TOTAL_INPUT_PATH)
if img is None:
    raise FileNotFoundError(f"Cannot load image: {TOTAL_INPUT_PATH}")

print("STEP: Original image")
imshow_rgb(img, title="Original Image (Full Page)", figsize=(12,8))

# 1) Convert to HSV
hsv = to_hsv(img)
print("STEP: HSV conversion (showing as HSV->RGB for display)")
imshow_rgb(hsv_to_display_rgb(hsv), title="HSV (display)")

# 2) Show S and V channels separately
print("STEP: Saturation and Value channels")
imshow_gray(hsv[:,:,1], title="Saturation (S) channel")
imshow_gray(hsv[:,:,2], title="Value (V) channel")

# 3) Threshold in HSV (S & V)
s_thresh = 30; v_thresh = 30
mask = mask_from_hsv(hsv, s_thresh, v_thresh)
print(f"STEP: HSV threshold mask (S>{s_thresh}, V>{v_thresh}) -> pixels= {np.sum(mask>0)}")
imshow_gray(mask, title=f"HSV Mask (S>{s_thresh}, V>{v_thresh})")

# 4) Morphological open (show)
m_open, m_close = morphological_clean(mask, open_k=3, close_k=5)
print("STEP: Morphological OPEN result")
imshow_gray(m_open, title="After Morphological Open")

# 5) Morphological close (show)
print("STEP: Morphological CLOSE result")
imshow_gray(m_close, title="After Morphological Close")

# 6) Contour extraction on closed mask
contours = contours_from_mask(m_close)
print(f"STEP: Contours extracted: {len(contours)}")
contours_overlay(img, contours, title="All Contours (from HSV mask)")

# 7) Filter contours -> circular blobs
blobs = contours_to_circles(contours, area_min=MIN_BLOB_AREA, circ_min=0.32)
print(f"STEP: Filtered circular blobs (area>={MIN_BLOB_AREA}, circ>=0.32): {len(blobs)}")
overlay_circles(img, blobs, title="Circular blobs (HSV + contour + circularity)")

# 8) Hough fallback (show)
hough_blobs = hough_circles_fallback(img, dp=1.2, minDist=28, param1=80, param2=26, minR=12, maxR=60)
print(f"STEP: Hough circles fallback candidates: {len(hough_blobs)}")
# overlay hough on original
overlay_circles(img, hough_blobs, title="Hough Circles (Fallback)")

# 9) Final combined view (prefer contour-based blobs, add Hough if needed)
final_blobs = blobs.copy()
# If contour-based found < expected and Hough has candidates, try to merge some non-duplicate
if len(final_blobs) < EXPECTED_COLS and len(hough_blobs) > 0:
    for hb in hough_blobs:
        hx,hy,hr,_,_ = hb
        too_close = False
        for b in final_blobs:
            if np.hypot(b[0]-hx, b[1]-hy) < 0.6*max(b[2], hr):
                too_close = True; break
        if not too_close:
            final_blobs.append(hb)
final_blobs = sorted(final_blobs, key=lambda b: b[0])
print(f"STEP: Final blobs used (combined): {len(final_blobs)}")
overlay_circles(img, final_blobs, title="Final Detected Wells (Full Image)")

STEP: Original image

STEP: HSV conversion (showing as HSV->RGB for display)

STEP: Saturation and Value channels

STEP: HSV threshold mask (S>30, V>30) -> pixels= 190745

STEP: Morphological OPEN result

STEP: Morphological CLOSE result

STEP: Contours extracted: 56

STEP: Filtered circular blobs (area>=60, circ>=0.32): 41

STEP: Hough circles fallback candidates: 263

STEP: Final blobs used (combined): 41

array([[[249, 248, 250],
        [249, 248, 250],
        [249, 248, 250],
        ...,
        [251, 249, 249],
        [251, 249, 249],
        [250, 248, 248]],

       [[249, 248, 250],
        [249, 248, 250],
        [249, 248, 250],
        ...,
        [254, 252, 252],
        [254, 252, 252],
        [254, 252, 252]],

       [[249, 248, 250],
        [249, 248, 250],
        [249, 248, 250],
        ...,
        [255, 254, 254],
        [255, 255, 255],
        [255, 255, 255]],

       ...,

       [[199, 194, 196],
        [199, 194, 196],
        [199, 194, 196],
        ...,
        [248, 250, 251],
        [248, 250, 251],
        [248, 250, 251]],

       [[199, 194, 196],
        [199, 194, 196],
        [199, 194, 196],
        ...,
        [248, 250, 251],
        [248, 250, 251],
        [248, 250, 251]],

       [[199, 194, 196],
        [199, 194, 196],
        [199, 194, 196],
        ...,
        [248, 250, 251],
        [248, 250, 251],
        [248, 250, 251]]], dtype=uint8)

# Cell 4: Cluster detected blobs into rows and create row crops

rows = cluster_rows_by_y_from_blobs(final_blobs)
print(f"Rows found by clustering: {len(rows)}")

# Show center lines for each row on the full image
overlay = img.copy()
for ridx, row_blobs in enumerate(rows):
    avg_y = int(np.mean([b[1] for b in row_blobs]))
    cv2.line(overlay, (0, avg_y), (overlay.shape[1], avg_y), (255,0,255), 2)
    # mark row centers
    for b in row_blobs:
        cv2.circle(overlay, (b[0], b[1]), 3, (0,255,255), -1)
imshow_rgb(overlay, title="Clustered Rows (center lines & blob centers)", figsize=(12,6))

# Create and show crops for each row (with margin)
row_crops = []
for i,row_blobs in enumerate(rows, start=1):
    xs = [b[0] for b in row_blobs]; ys = [b[1] for b in row_blobs]; rs = [b[2] for b in row_blobs]
    if len(xs) == 0:
        continue
    x_min = max(0, int(min(xs) - max(rs) - 20))
    x_max = min(img.shape[1], int(max(xs) + max(rs) + 20))
    y_min = max(0, int(min(ys) - max(rs) - 20))
    y_max = min(img.shape[0], int(max(ys) + max(rs) + 20))
    crop = img[y_min:y_max, x_min:x_max].copy()
    row_crops.append((i, crop, (x_min,y_min,x_max,y_max)))
    imshow_rgb(crop, title=f"Row crop {i} (blobs: {len(row_blobs)})", figsize=(10,3))

print(f"Created {len(row_crops)} row crops.")

Rows found by clustering: 3

Created 3 row crops.

# Cell 5: For each row crop, run stepwise detection and extract RGB + S_mean
all_trial_dfs = []

for (rid, crop, bbox) in row_crops:
    print(f"\n=== Processing row crop {rid} ===")
    imshow_rgb(crop, title=f"Row {rid} - Original crop", figsize=(10,3))

    hsv_c = to_hsv(crop)
    imshow_rgb(hsv_to_display_rgb(hsv_c), title=f"Row {rid} - HSV (display)", figsize=(10,3))
    imshow_gray(hsv_c[:,:,1], title=f"Row {rid} - S channel")
    imshow_gray(hsv_c[:,:,2], title=f"Row {rid} - V channel")

    # threshold
    mask_c = mask_from_hsv(hsv_c, s_thresh=30, v_thresh=30)
    imshow_gray(mask_c, title=f"Row {rid} - HSV mask (S>30,V>30)")

    # morphological open & close
    m_open_c, m_close_c = morphological_clean(mask_c, open_k=3, close_k=5)
    imshow_gray(m_open_c, title=f"Row {rid} - Morph open")
    imshow_gray(m_close_c, title=f"Row {rid} - Morph close")

    # contours
    cnts_c = contours_from_mask(m_close_c)
    print(f"Row {rid}: contours found = {len(cnts_c)}")
    contours_overlay(crop, cnts_c, title=f"Row {rid} - All contours")

    # circularity filtering
    blobs_c = contours_to_circles(cnts_c, area_min=MIN_BLOB_AREA, circ_min=0.3)
    print(f"Row {rid}: blobs after area/circularity filter = {len(blobs_c)}")
    overlay_circles(crop, blobs_c, title=f"Row {rid} - Circular blobs (filtered)")

    # fallback: add hough if needed
    if len(blobs_c) < EXPECTED_COLS:
        hough_c = hough_circles_fallback(crop, dp=1.2, minDist=24, param1=80, param2=26, minR=12, maxR=60)
        print(f"Row {rid}: Hough fallback candidates = {len(hough_c)}")
        overlay_circles(crop, hough_c, title=f"Row {rid} - Hough fallback")
        # merge non-duplicates
        for hb in hough_c:
            hx,hy,hr,_,_ = hb
            too_close = False
            for b in blobs_c:
                if np.hypot(b[0]-hx, b[1]-hy) < 0.6*max(b[2], hr):
                    too_close = True; break
            if not too_close:
                blobs_c.append(hb)
        blobs_c = sorted(blobs_c, key=lambda b: b[0])
        overlay_circles(crop, blobs_c, title=f"Row {rid} - Combined blobs")

    # final selection: at most EXPECTED_COLS, take largest if more
    if len(blobs_c) > EXPECTED_COLS:
        blobs_c = sorted(blobs_c, key=lambda b: b[3], reverse=True)[:EXPECTED_COLS]
        blobs_c = sorted(blobs_c, key=lambda b: b[0])

    overlay_circles(crop, blobs_c, title=f"Row {rid} - Final blobs (annotated)")

    # If number of blobs is not EXPECTED_COLS, warn (we will still extract what we have)
    if len(blobs_c) != EXPECTED_COLS:
        print(f"Warning: row {rid} has {len(blobs_c)} detected wells (expected {EXPECTED_COLS}).")

    # Extract inner RGB + S for each detected blob (left->right)
    records = []
    for idx, b in enumerate(blobs_c[:EXPECTED_COLS]):   # use only first EXPECTED_COLS if more
        x,y,r,area,circ = b
        Rval,Gval,Bval,RGBmean,Smean = extract_inner_rgb_and_saturation(crop, hsv_c, x, y, r)
        # map concentrations skipping control at index 0
        if idx == 0:
            conc = None
        else:
            conc = CONCENTRATIONS[idx-1] if idx-1 < len(CONCENTRATIONS) else None
        records.append({
            "Trial": rid,
            "Well": idx+1,
            "Concentration": conc,
            "R": Rval, "G": Gval, "B": Bval, "RGB_mean": RGBmean, "S_mean": Smean
        })
    df_row = pd.DataFrame(records, columns=["Trial","Well","Concentration","R","G","B","RGB_mean","S_mean"])
    print(f"Row {rid} extracted dataframe:")
    display(df_row)
    all_trial_dfs.append((rid, df_row))

=== Processing row crop 1 ===

Row 1: contours found = 17

Row 1: blobs after area/circularity filter = 13

Row 1 extracted dataframe:

	Trial	Well	Concentration	R	G	B	RGB_mean	S_mean
0	1	1	NaN	241.173975	229.983027	183.999057	218.385353	60.507779
1	1	2	0.5	178.925067	215.213769	200.073792	198.070876	43.051350
2	1	3	1.0	168.120198	210.516546	200.736402	193.124382	51.408520
3	1	4	2.0	153.456930	205.127614	202.250266	186.944937	64.357320
4	1	5	3.0	152.904998	203.729529	199.699752	185.444760	63.658986
5	1	6	4.0	145.106631	201.277574	202.266578	182.883594	72.790736
6	1	7	5.0	142.227224	200.953563	203.613258	182.264682	77.064516
7	1	8	6.0	137.264445	199.589507	205.485289	180.779747	84.672102
8	1	9	7.0	133.848384	199.744085	208.281573	180.624681	91.112296
9	1	10	8.0	129.863446	196.822366	204.875998	177.187270	93.375428
10	1	11	9.0	133.673442	200.028657	210.800733	181.500944	93.210263
11	1	12	10.0	122.045728	196.789791	210.133641	176.323053	106.878766

=== Processing row crop 2 ===

Row 2: contours found = 30

Row 2: blobs after area/circularity filter = 19

Row 2 extracted dataframe:

	Trial	Well	Concentration	R	G	B	RGB_mean	S_mean
0	2	1	NaN	239.892946	227.476072	178.818859	215.395959	64.935838
1	2	2	0.5	180.392769	215.974123	201.439915	199.268935	42.068061
2	2	3	1.0	170.055915	211.914797	202.179536	194.716749	50.430202
3	2	4	2.0	156.432471	206.071606	202.734846	188.412974	61.520028
4	2	5	3.0	152.571476	204.901033	204.483172	187.318560	65.799067
5	2	6	4.0	146.088501	202.211592	205.005921	184.435338	73.434403
6	2	7	5.0	142.261720	200.692214	203.541378	182.165104	76.878516
7	2	8	6.0	140.119437	200.592621	206.357931	182.356663	81.903766
8	2	9	7.0	132.351549	200.019993	212.301233	181.557592	96.036654
9	2	10	8.0	133.995392	199.863169	207.197802	180.352121	90.180078
10	2	11	9.0	132.396348	199.330544	208.765310	180.164067	93.282617
11	2	12	10.0	123.890464	196.753633	210.687345	177.110481	105.067706

=== Processing row crop 3 ===

Row 3: contours found = 28

Row 3: blobs after area/circularity filter = 21

Row 3 extracted dataframe:

	Trial	Well	Concentration	R	G	B	RGB_mean	S_mean
0	3	1	NaN	239.807870	227.497341	180.710032	216.005081	62.852180
1	3	2	0.5	180.209966	216.487258	204.551541	200.416255	42.793838
2	3	3	1.0	168.604750	211.494505	203.212691	194.437315	51.770649
3	3	4	2.0	159.248052	207.575569	204.976317	190.599979	59.493300
4	3	5	3.0	155.056315	205.210930	203.129957	187.799067	62.549817
5	3	6	4.0	150.979794	203.783410	205.220844	186.661349	67.695498
6	3	7	5.0	143.573910	202.886211	209.749380	185.403167	80.496987
7	3	8	6.0	139.904365	200.138954	206.614129	182.219149	82.473176
8	3	9	7.0	137.621411	201.537398	209.950727	183.036512	87.902871
9	3	10	8.0	133.723502	199.121943	209.900035	180.915160	92.558667
10	3	11	9.0	131.946118	200.461538	211.651188	181.352948	96.031549
11	3	12	10.0	130.914215	198.462602	209.472173	179.616330	95.658986

# Cell 6: Collect valid trials (exactly EXPECTED_COLS and S_mean >= threshold for all wells except control)
valid_trials = []
for (rid, df_row) in all_trial_dfs:
    if len(df_row) < EXPECTED_COLS:
        print(f"Row {rid} skipped: only {len(df_row)} wells detected.")
        continue
    # check S_mean for wells (skip control idx 0)
    svals = df_row["S_mean"].values
    # consider all wells present and non-NaN
    if np.any(np.isnan(svals)):
        print(f"Row {rid} skipped due to NaN S values.")
        continue
    # verify S >= threshold for all non-control wells
    if np.all(svals >= SATURATION_THRESHOLD):
        valid_trials.append((rid, df_row))
    else:
        print(f"Row {rid} skipped due to low saturation in some wells (S mean min={np.min(svals):.2f}).")

print(f"\nValid trials count: {len(valid_trials)}")
combined_rows = []
for (tid, df_row) in valid_trials:
    # drop control (Well==1) when building concentration mapping
    df_use = df_row[df_row["Well"] != 1].copy().reset_index(drop=True)
    # ensure concentrations column is filled (should be from previous step)
    combined_rows.append((tid, df_use))

if len(combined_rows) == 0:
    print("No valid trials to analyze further.")
else:
    for tid, df_use in combined_rows:
        print(f"\n--- Trial {tid} (used for calibration) ---")
        display(df_use)

Valid trials count: 3

--- Trial 1 (used for calibration) ---

	Trial	Well	Concentration	R	G	B	RGB_mean	S_mean
0	1	2	0.5	178.925067	215.213769	200.073792	198.070876	43.051350
1	1	3	1.0	168.120198	210.516546	200.736402	193.124382	51.408520
2	1	4	2.0	153.456930	205.127614	202.250266	186.944937	64.357320
3	1	5	3.0	152.904998	203.729529	199.699752	185.444760	63.658986
4	1	6	4.0	145.106631	201.277574	202.266578	182.883594	72.790736
5	1	7	5.0	142.227224	200.953563	203.613258	182.264682	77.064516
6	1	8	6.0	137.264445	199.589507	205.485289	180.779747	84.672102
7	1	9	7.0	133.848384	199.744085	208.281573	180.624681	91.112296
8	1	10	8.0	129.863446	196.822366	204.875998	177.187270	93.375428
9	1	11	9.0	133.673442	200.028657	210.800733	181.500944	93.210263
10	1	12	10.0	122.045728	196.789791	210.133641	176.323053	106.878766

--- Trial 2 (used for calibration) ---

	Trial	Well	Concentration	R	G	B	RGB_mean	S_mean
0	2	2	0.5	180.392769	215.974123	201.439915	199.268935	42.068061
1	2	3	1.0	170.055915	211.914797	202.179536	194.716749	50.430202
2	2	4	2.0	156.432471	206.071606	202.734846	188.412974	61.520028
3	2	5	3.0	152.571476	204.901033	204.483172	187.318560	65.799067
4	2	6	4.0	146.088501	202.211592	205.005921	184.435338	73.434403
5	2	7	5.0	142.261720	200.692214	203.541378	182.165104	76.878516
6	2	8	6.0	140.119437	200.592621	206.357931	182.356663	81.903766
7	2	9	7.0	132.351549	200.019993	212.301233	181.557592	96.036654
8	2	10	8.0	133.995392	199.863169	207.197802	180.352121	90.180078
9	2	11	9.0	132.396348	199.330544	208.765310	180.164067	93.282617
10	2	12	10.0	123.890464	196.753633	210.687345	177.110481	105.067706

--- Trial 3 (used for calibration) ---

	Trial	Well	Concentration	R	G	B	RGB_mean	S_mean
0	3	2	0.5	180.209966	216.487258	204.551541	200.416255	42.793838
1	3	3	1.0	168.604750	211.494505	203.212691	194.437315	51.770649
2	3	4	2.0	159.248052	207.575569	204.976317	190.599979	59.493300
3	3	5	3.0	155.056315	205.210930	203.129957	187.799067	62.549817
4	3	6	4.0	150.979794	203.783410	205.220844	186.661349	67.695498
5	3	7	5.0	143.573910	202.886211	209.749380	185.403167	80.496987
6	3	8	6.0	139.904365	200.138954	206.614129	182.219149	82.473176
7	3	9	7.0	137.621411	201.537398	209.950727	183.036512	87.902871
8	3	10	8.0	133.723502	199.121943	209.900035	180.915160	92.558667
9	3	11	9.0	131.946118	200.461538	211.651188	181.352948	96.031549
10	3	12	10.0	130.914215	198.462602	209.472173	179.616330	95.658986

# Cell 7: Per-trial R-only analysis and plotting

def fit_inverse_model_and_report(df_trial, trial_id, deg_forward=2, deg_inverse=2):
    # df_trial should have Concentration and R columns, no control
    C = df_trial["Concentration"].values
    R = df_trial["R"].values.reshape(-1,1)

    # forward model R = f(C) for plotting
    poly_f = PolynomialFeatures(deg_forward, include_bias=False)
    Xf = poly_f.fit_transform(C.reshape(-1,1))
    model_f = LinearRegression().fit(Xf, R.ravel())

    # inverse model C = g(R) for prediction
    poly_i = PolynomialFeatures(deg_inverse, include_bias=False)
    Xi = poly_i.fit_transform(R)
    model_i = LinearRegression().fit(Xi, C)

    # predictions & metrics
    C_pred_from_R = model_i.predict(poly_i.transform(R))
    r2 = r2_score(C, C_pred_from_R)
    mae = mean_absolute_error(C, C_pred_from_R)
    rmse = math.sqrt(mean_squared_error(C, C_pred_from_R))

    # plot R vs C and forward fit
    xs = np.linspace(min(C), max(C), 200).reshape(-1,1)
    ys = model_f.predict(poly_f.transform(xs))

    plt.figure(figsize=(6,4))
    plt.scatter(C, R, label="observed (R)", color='blue')
    plt.plot(xs, ys, 'r-', label="forward fit R=f(C)")
    plt.xlabel("Concentration")
    plt.ylabel("R channel")
    plt.title(f"Trial {trial_id}: R vs Concentration (deg_fwd={deg_forward})")
    plt.legend(); plt.grid(True)
    plt.show()

    # Print metrics and per-well predictions
    print(f"Trial {trial_id} metrics: R2={r2:.4f}, MAE={mae:.4f}, RMSE={rmse:.4f}")
    df_pred = pd.DataFrame({
        "Well": df_trial["Well"].values,
        "True_Conc": C,
        "R_value": R.ravel(),
        "Pred_Conc_from_R": C_pred_from_R
    })
    display(df_pred)
    return {"r2": r2, "mae": mae, "rmse": rmse, "model_inv": model_i, "poly_inv": poly_i, "model_fwd": model_f, "poly_fwd": poly_f}

# Run for each valid trial
results = {}
for tid, df_use in combined_rows:
    results[tid] = fit_inverse_model_and_report(df_use, tid, deg_forward=2, deg_inverse=2)

Trial 1 metrics: R2=0.9541, MAE=0.4072, RMSE=0.6629

	Well	True_Conc	R_value	Pred_Conc_from_R
0	2	0.5	178.925067	0.336259
1	3	1.0	168.120198	0.944719
2	4	2.0	153.456930	2.784838
3	5	3.0	152.904998	2.876912
4	6	4.0	145.106631	4.354717
5	7	5.0	142.227224	4.983883
6	8	6.0	137.264445	6.173993
7	9	7.0	133.848384	7.070935
8	10	8.0	129.863446	8.197355
9	11	9.0	133.673442	7.118575
10	12	10.0	122.045728	10.657813

Trial 2 metrics: R2=0.9693, MAE=0.4232, RMSE=0.5416

	Well	True_Conc	R_value	Pred_Conc_from_R
0	2	0.5	180.392769	0.434614
1	3	1.0	170.055915	0.862787
2	4	2.0	156.432471	2.401778
3	5	3.0	152.571476	3.039515
4	6	4.0	146.088501	4.310573
5	7	5.0	142.261720	5.178660
6	8	6.0	140.119437	5.702810
7	9	7.0	132.351549	7.833232
8	10	8.0	133.995392	7.352332
9	11	9.0	132.396348	7.819912
10	12	10.0	123.890464	10.563786

Trial 3 metrics: R2=0.9898, MAE=0.2424, RMSE=0.3131

	Well	True_Conc	R_value	Pred_Conc_from_R
0	2	0.5	180.209966	0.580173
1	3	1.0	168.604750	0.986845
2	4	2.0	159.248052	2.059969
3	5	3.0	155.056315	2.756510
4	6	4.0	150.979794	3.561975
5	7	5.0	143.573910	5.348393
6	8	6.0	139.904365	6.387970
7	9	7.0	137.621411	7.086363
8	10	8.0	133.723502	8.370341
9	11	9.0	131.946118	8.994142
10	12	10.0	130.914215	9.367319